#include <iostream>
#include <fstream>
#include <string>
#include <math.h>
#include <time.h>
#include <float.h>
#include "SDL.h"
#include "SDL_opengl.h"
#include "Vector3D.h"
#include "Quaternion.h"
#include "metatable.h"
const int width=800;
const int height=600;
bool keys[SDLK_LAST]={0};
//#define DEBUG
#define NBALLS (5)
#define NCUBES (20)
#define CUBE_SZ (0.2)
int MIN_THRES=0;
#define SCENE_SZ (NCUBES*CUBE_SZ)
#ifdef M_PI
#undef M_PI
#endif
#define M_PI (3.14159265358979323846)
Vector3D light_pos(0,10,0);
class MetaSurface{
protected:
	int n;
	Vector3D* x;
public:
	MetaSurface(){	}
	MetaSurface(int n){
		this->n=n;
		this->x=new Vector3D[n];
	}
	~MetaSurface(){
		delete[] x;
	}
	virtual const double potential(Vector3D& p) const{
		if(p[0]>SCENE_SZ || p[1]>SCENE_SZ || p[2]>SCENE_SZ || p[0]<0 || p[1]<0 || p[2]<0)
			return -1;
		double sum=0;
		for(int i=0;i<n;i++){
			const double _r = (p-x[i]).squareNorm();
			sum+=1.0/(1.0+_r);
		}
		return sum;
	}
	virtual void update(){ }
	const bool outside(Vector3D p) const{ return potential(p) < MIN_THRES; }
	Vector3D intersection(Vector3D& p1, Vector3D& p2) const{
		const double v1=potential(p1);
		const double denom=potential(p2)-v1;
		if(fabs(denom)<1e-5) return p1;	//Prevent divide by zero
		const double m=(MIN_THRES-v1)/denom;
		return p1+m*(p2-p1);			//Linearly interpolate between the points
	}
};
class MetaFile : public MetaSurface{
private:
	double* k;
	int kMAX;
public:
	MetaFile(int p, istream& in) : MetaSurface(p){
		this->k=new double[this->n];
		int lats=(int)sqrt(n);
		int longs=n/lats;
		kMAX=0;
		for(int i=0;i<this->n;i++){
			double tmp=k[i]=0;
			for(int j=0;j<this->n;j++){
				in>>tmp;
				k[i]+=tmp*tmp;
			}
			k[i]=sqrt(k[i]);
			if(k[i]>kMAX) kMAX=k[i];
			cout<<"k["<<i<<"]="<<k[i]<<endl;
		}
		for(int i=0;i<lats;i++){	//Start at the bottom and iterate through the latitudes
			const double z1=sin(2*M_PI*i/lats);	//Figure out the Z-coord and radius of the two lats
			const double zr1=cos(2*M_PI*i/lats);
			for(int j=0;j<longs;j++){
				const double X=sin(2*M_PI*j/longs);	//Find the X,Y point of the circles
				const double Y=cos(2*M_PI*j/longs);
				x[i*lats+j]=Vector3D(X*zr1,Y*zr1,z1)*(k[i*lats+j]/kMAX)+Vector3D(SCENE_SZ/2.0,SCENE_SZ/2.0,SCENE_SZ/2.0);		//and scale them according to the layer radius, and draw those points
			}
		}
		cout<<SCENE_SZ/2.0<<endl;
		for(int i=0;i<n;i++) cout<<x[i]<<endl;
	}
	const double potential(Vector3D& p) const{
		if(p[0]>SCENE_SZ || p[1]>SCENE_SZ || p[2]>SCENE_SZ || p[0]<0 || p[1]<0 || p[2]<0)
			return -1;
		double sum=0;
		for(int i=0;i<n;i++){
			const double _r = (p-x[i]).squareNorm();
			sum+=k[i]/(_r);
		}
		return sum;
	}
	~MetaFile(){
		delete[] k;
	}
};
class MarchCubes{
private:
	MetaSurface* s;
	bool* visted;
	const int getIndex(Vector3D* p) const{	//Edge intersection table index.
		int x=0;
		for(int i=0;i<8;i++)
			if(s->outside(p[i]))
				x|=(1<<i);
		return x;
	}
	void triangulate(Vector3D* c, const int index){	//Generates the triangles for a given voxel
		Vector3D verts[12];
		const unsigned int edge=edgeTable[index];
		if(!edge) return;
		if(edge & 1   ) verts[ 0] = s->intersection(c[0], c[1]);
		if(edge & 2   ) verts[ 1] = s->intersection(c[1], c[2]);
		if(edge & 4   ) verts[ 2] = s->intersection(c[2], c[3]);
		if(edge & 8   ) verts[ 3] = s->intersection(c[3], c[0]);
		if(edge & 16  ) verts[ 4] = s->intersection(c[4], c[5]);
		if(edge & 32  ) verts[ 5] = s->intersection(c[5], c[6]);
		if(edge & 64  ) verts[ 6] = s->intersection(c[6], c[7]);
		if(edge & 128 ) verts[ 7] = s->intersection(c[7], c[4]);
		if(edge & 256 ) verts[ 8] = s->intersection(c[0], c[4]);
		if(edge & 512 ) verts[ 9] = s->intersection(c[1], c[5]);
		if(edge & 1024) verts[10] = s->intersection(c[2], c[6]);
		if(edge & 2048) verts[11] = s->intersection(c[3], c[7]);
		for(int i=0; triTable[index][i]!=-1;i+=3, nVerts+=3){
			mesh[ nVerts ] = verts[triTable[index][i+2]];
			mesh[nVerts+1] = verts[triTable[index][i+1]];
			mesh[nVerts+2] = verts[triTable[index][ i ]];
			Vector3D n = (mesh[nVerts+1]-mesh[nVerts]).crossProduct(mesh[nVerts+2]-mesh[nVerts]);
			n.normalize();
			normal[ nVerts ] = n;
			normal[nVerts+1] = n;
			normal[nVerts+2] = n;
		}
	}
public:
	Vector3D* mesh;
	Vector3D* normal;
	int nVerts;
	MarchCubes(MetaSurface* s){
		this->s=s;
		mesh=new Vector3D[15*(int)(NCUBES*NCUBES*NCUBES)];
		normal=new Vector3D[15*(int)(NCUBES*NCUBES*NCUBES)];
		visted=new bool[(int)(NCUBES*NCUBES*NCUBES)];
		reset();
	}
	~MarchCubes(){
		delete[] mesh;
		delete[] visted;
		s=NULL;
	}
	void reset(){
		memset(visted,0,(int)(NCUBES*NCUBES*NCUBES)*sizeof(bool));
		nVerts=0;
	}
	void generate(const double x, const double y, const double z){
		if(x>SCENE_SZ || y>SCENE_SZ || z>SCENE_SZ || x<0 || y<0 || z<0) return; //Out of bounds.
		const int index = ((z) /CUBE_SZ) * NCUBES * NCUBES + ((y) / CUBE_SZ) * NCUBES + ((x) / CUBE_SZ);
		if(visted[index]) return;
		#ifdef DEBUG
		cout<<"COMPUTING CUBE "<<index<<endl;
		#endif
		visted[index]=true;
		Vector3D cube[8];
		cube[0]=Vector3D(x,y,z+CUBE_SZ);
		cube[1]=Vector3D(x+CUBE_SZ,y,z+CUBE_SZ);
		cube[2]=Vector3D(x+CUBE_SZ,y,z);
		cube[3]=Vector3D(x,y,z);
		cube[4]=Vector3D(x,y+CUBE_SZ,z+CUBE_SZ);
		cube[5]=Vector3D(x+CUBE_SZ,y+CUBE_SZ,z+CUBE_SZ);
		cube[6]=Vector3D(x+CUBE_SZ,y+CUBE_SZ,z);
		cube[7]=Vector3D(x,y+CUBE_SZ,z);
		const int ind = getIndex(cube);
		#ifdef DEBUG
		cout<<"ind:"<<ind<<endl;
		#endif
		//if(!edgeTable[ind]) return;	//For optimization, should be implemented
		triangulate(cube, ind);
		//NEED TO MAKE THIS ITERATIVE
		generate(x-CUBE_SZ, y, z);
		generate(x+CUBE_SZ, y, z);
		generate(x, y-CUBE_SZ, z);
		generate(x, y+CUBE_SZ, z);
		generate(x, y, z-CUBE_SZ);
		generate(x, y, z+CUBE_SZ);
	}
};
MetaFile* s;
MarchCubes* mc;
inline void DrawScene(){
	static double angle=0.0;
	gluLookAt(0,0,-10,0,0,0,0,1,0);
	glPushMatrix();
	glRotated(angle,0,1,0);
	glTranslated(-SCENE_SZ/2.0,-SCENE_SZ/2.0,-SCENE_SZ/2.0);
	angle+=0.5;
	glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
	glEnable(GL_CULL_FACE);
	glBegin(GL_TRIANGLES);
	for(int i=0;i<mc->nVerts;i+=3){
		glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, (GLfloat[4]){mc->normal[i][0],mc->normal[i][1],mc->normal[i][2],0.75});
		glNormal3d((mc->normal[i+2])[0],(mc->normal[i+2])[1],(mc->normal[i+2])[2]);
		glVertex3d((mc->mesh[i+2])[0],(mc->mesh[i+2])[1],(mc->mesh[i+2])[2]);
		glNormal3d((mc->normal[i+1])[0],(mc->normal[i+1])[1],(mc->normal[i+1])[2]);
		glVertex3d((mc->mesh[i+1])[0],(mc->mesh[i+1])[1],(mc->mesh[i+1])[2]);
		glNormal3d((mc->normal[i])[0],(mc->normal[i])[1],(mc->normal[i])[2]);
		glVertex3d((mc->mesh[i])[0],(mc->mesh[i])[1],(mc->mesh[i])[2]);
	}
	glEnd();
	glPushMatrix();
	glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
	glDisable(GL_CULL_FACE);
	glScaled(SCENE_SZ,SCENE_SZ,SCENE_SZ);
	glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, (GLfloat[4]){0,0,0.5,1});
	glBegin(GL_QUADS);
		glVertex3d(0,0,0);	//Left side
		glVertex3d(0,0,1);
		glVertex3d(0,1,1);
		glVertex3d(0,1,0);

		glVertex3d(0,0,0);	//Front
		glVertex3d(0,1,0);
		glVertex3d(1,1,0);
		glVertex3d(1,0,0);

		glVertex3d(0,0,0);	//Bottom
		glVertex3d(1,0,0);
		glVertex3d(1,0,1);
		glVertex3d(0,0,1);

		glVertex3d(1,0,0);	//Right side
		glVertex3d(1,1,0);
		glVertex3d(1,1,1);
		glVertex3d(1,0,1);

		glVertex3d(0,0,1);	//Back
		glVertex3d(1,0,1);
		glVertex3d(1,1,1);
		glVertex3d(0,1,1);

		glVertex3d(0,1,0);	//Top
		glVertex3d(0,1,1);
		glVertex3d(1,1,1);
		glVertex3d(1,1,0);
	glEnd();
	glPopMatrix();
}
int main (int argc, char ** argv) {
	if(argc<3){ printf("Usage: %s n filename\n",argv[0]); return 1; }
	else if((MIN_THRES=atoi(argv[1]))< 0){ printf("Usage: %s n\n",argv[0]); return 1; }
	srand((unsigned long)time(NULL));			//Setup random
	if (SDL_Init(SDL_INIT_VIDEO) < 0) exit(1);		//Init SDL
	SDL_WM_SetCaption("Metaballs3D", "Metaballs3D");			//Title and icon(?) captions
	const SDL_VideoInfo *video = SDL_GetVideoInfo();	//Get current video info (for pixel format later)
	if (!video) { SDL_Quit(); exit(1); }
	SDL_Surface *screen = SDL_SetVideoMode(width, height, video->vfmt->BitsPerPixel, SDL_OPENGL);	//Setup video
	if (!screen) { SDL_Quit(); exit(1); }
	glViewport(0, 0, width, height);			//Setup the entire window as a viewport
	
	glDisable(GL_CULL_FACE);
	//glCullFace(GL_BACK);
	glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);
	glDepthFunc(GL_LEQUAL);
	//glEnable(GL_DEPTH_TEST);
	glDisable(GL_DEPTH_TEST);  //Might want to for various reasons (2D, alpha etc).
	glDepthMask(GL_TRUE);
	glClearDepth(1.0);
	glClearColor(0.0, 0.0, 0.0, 0.0);
	glShadeModel(GL_SMOOTH);

	//Blending!
	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);

	
	//Lighting
	glLightfv(GL_LIGHT1, GL_AMBIENT, (GLfloat[4]){0.5,0.5,0.5,1.0});
	glLightfv(GL_LIGHT1, GL_DIFFUSE, (GLfloat[4]){1.0,1.0,1.0,1.0});
	glLightfv(GL_LIGHT1, GL_POSITION,(GLfloat[4]){light_pos[0],light_pos[1],light_pos[2],1.0});
	glEnable(GL_LIGHT1);
	glEnable(GL_LIGHTING);
	glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, (GLfloat[4]){1,0,0,0.5});
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(45.0f,(GLfloat)width/(GLfloat)height,1.0f,100.0f); //3D
	//glOrtho(0,width,height,0,0,1);	//2D
	glMatrixMode (GL_MODELVIEW);
	glLoadIdentity();
	SDL_GL_SetAttribute(SDL_GL_RED_SIZE, 5);   		// Set Pixel Format
	SDL_GL_SetAttribute(SDL_GL_GREEN_SIZE, 5);
	SDL_GL_SetAttribute(SDL_GL_BLUE_SIZE, 5);
	SDL_GL_SetAttribute(SDL_GL_DEPTH_SIZE, 16);
	SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1);
	//Other init
	ifstream in(argv[2]);
	int n;
	in>>n;
	s=new MetaFile(n,in);
	mc=new MarchCubes(s);
	bool done=false;
	while(!done){
		SDL_Event event;
		if(SDL_PollEvent(&event)){	//Get an event
			switch (event.type){
				case SDL_KEYDOWN:
					switch(event.key.keysym.sym){
						case SDLK_ESCAPE:
							done=true; continue;
						default:	//Any keys pressed, it'll be true in the global array.
							keys[event.key.keysym.sym]=true; break;
					}
					break;
				case SDL_KEYUP:			//If let up, it'll be false in the array
					keys[event.key.keysym.sym]=false; break;
				case SDL_MOUSEBUTTONDOWN:
					break;
				case SDL_QUIT:			//Window quits
					done=true; continue;
				default:
					break;
			}
		}
		//Key check here.
		s->update();
		mc->reset();
		mc->generate(0,0,0);

		glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);	//Clear the screen
		glLoadIdentity();		//Reset model matrix.
		DrawScene();
		SDL_GL_SwapBuffers();		//Swap the buffer to the screen
	}
	delete s;
	delete mc;
	SDL_Quit();				//Kill and free all SDL engine information
	return 0;
}
